active_sites=[1 2 4 7 8];
passive_sites=[3 5 6 9];
control_sites=10;

mean_resp_by_site=nan(83,9);
aaa=[];
for curr_exp = 1:9
    % % time_point = 1;
    for time_point = 1:4
        cur_act=[];
        for layer = 1:4
            cur_act(end+1:end+size(proj_meta(curr_exp).rd(layer,time_point).act(:,1:7500),1),:) = proj_meta(curr_exp).rd(layer,time_point).act(:,1:7500);
        end
        % generate binary activity during feedback.
        velM_bi = proj_meta(curr_exp).rd(1,time_point).velM_smoothed(1:7500);
        velM_bi = -velM_bi; % make the velocity of forward-running positive.
        velM_bi(velM_bi > 0.001) = 1;
        velM_bi(velM_bi <= 0.001) = 0;
        
        % generate the frames of ps onset.
        ps_ons = find(diff(proj_meta(curr_exp).rd(1,time_point).ps_id(1:7500))>1);
        % remove early & late ps.
        ps_ons(ps_ons > 7425) = [];
        ps_ons(ps_ons < 6) = [];
        
        rps_nbr = 0; % number of real ps
        rps_id = [];
        for psid = 1:length(ps_ons)
            if sum(velM_bi(ps_ons(psid)-5:ps_ons(psid)+15)) >= 15
                rps_nbr = rps_nbr+1;
                rps_id(rps_nbr) = psid;
            end
        end

        for ind = 1:length(rps_id)
            rps_frame_id = ps_ons(rps_id(ind));
            tmp = mean(cur_act(:,rps_frame_id-7:rps_frame_id+75),1);            
            aaa(curr_exp,time_point,ind)=mean(tmp(39:53))-mean(tmp(1:7));
        end        
    end
end

aaa(aaa==0)=nan;

figure; hold on
plot(smooth2(squeeze(nanmean(nanmean(aaa(active_sites,[1:4],:),1),2)),3))
plot(smooth2(squeeze(nanmean(nanmean(aaa(passive_sites,[1:4],:),1),2)),3),'r')